Continuous Time Quantum Monte Carlo Method for Fermions: 
Beyond Auxiliary Field Framework 
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Numerically exact continuous-time Quantum Monte Carlo algorithm for finite fermionic systems 
with non-local interactions is proposed. The scheme is particularly applicable for general multi- 
band time-dependent correlations since it does not invoke Hubbard- Stratonovich transformation. 
The present determinantal grand-canonical method is based on a stochastic series expansion for 
the partition function in the interaction representation. The results for the Green function and for 
the time-dependent susceptibility of multi-orbital super-symmetric impurity model with a spin-flip 
interaction are presented. 
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Quantum Monte Carlo (QMC) tools for fermionic sys- 
tems appeared more than 20 years ago 0,|3ilHB and are 
nowdays vital for a wide range of fields, like the physics of 
correlated materials, quantum chemistry and nanoelec- 
tronics. Although the first programs were developed for 
model Hamiltonians with local interaction, many-particle 
action of a very general form stays behind the real sys- 
tems. For example all matrix elements of the interaction 
do not vanish in the problems of quantum chemistry 5] 
and solid state physics 0]. Dynamical mean-field the- 
ory (DMFT) for correlated materials brings a non- 
trivial bath Green function on the scene, and its exten- 
sion deals with an interaction which is non-local in 
time. An off-diagonal exchange term can be responsible 
for the correlated superconductivity in doped fuUerens 
. It is worth to note in general that exchange is often 
of an indirect origin (like super-exchange) and the ex- 
change terms are therefore retarded. New developments 
|lOj| clearly urge an invention of essentially different type 
of QMC scheme suitable for non-local, time-dependent 
interaction. 

The determinantal grand-canonical auxiliary-field 
scheme Q, |^ IE Q is commonly used for the inter- 
acting fermions, because other known QMC schemes 
(like stochastic series expansion in powers of Hamil- 
tonian or worm algorithms suffer an unac- 

ceptably bad sign problem for this case. Two points 
are essential for the approach: first, the imaginary 
time is artificially discretized, and then the Hubbard- 
Stratonovich transformation |l3| is performed to decou- 
ple the fermionic degrees of freedom. After the decou- 
pling, fermions can be integrated out, and Monte Carlo 
sampling should be performed in the s_pace of auxil- 
iary Hubbard-Stratonovich fields. Hirsh [3] proposed to 
use discrete Hubbard-Stratonovich transformation to im- 
prove the original scheme; this is now a standard method 
for simulations of lattice and impurity quantum prob- 
lems. For relatively small clusters, and in particular for 



DMFT, the sign problem is not crucial in this method 
il4]. The number of auxiliary field is linear (quadratic) 
in the number of atoms for the case of local (nonlocal) 
interaction. 

The time discretization leads in a systematic error 
of the result. For for bosonic quantum systems, con- 
tinuous time loop algorithm 1 151 . worm diagrammatic 
world line Monte Carlo scheme \z\ and continuous time 
path-integral QMC 16] overcame this issue. Recently a 
continuous- time modification of the fermionic QMC al- 
gorithm was proposed 1^. It is based on a series ex- 
pansion for the partition function in the powers of inter- 
action. The scheme is free of time-discretization errors, 
but the Hubbard-Stratonovich transformation is still in- 
voked. Therefore the number of auxiliary fields scales 
similarly to the discrete scheme. This scheme is devel- 
oped for local interaction only. 

Besides the time-discretization problem, the non- 
locality of interaction hampers the calculation in the ex- 
isting schemes, because it is hard to simulate systems 
with a large number of auxiliary spins. Further, the dis- 
crete Hubbard-Stratonovich transformation is not suit- 
able for non-local in time interactions. One needs to use 
continuous dispersive bosonic fields |^ for this case, that 
makes the simulation even harder. 

In this Letter we present a novel numerically ex- 
act continuous-time fermionic QMC algorithm. This is 
the first QMC scheme that do not invoke any type of 
Hubbard-Stratonovich transformations and therefore op- 
erates natively with non-local in space and time inter- 
actions. The scheme is free of systematic errors due 
to direct operations with continuous time expansion of 
the partition function. Numerical results for a super- 
symmetric two band impurity model with spin-flip, time 
dependent non-local interactions show an advantage and 
a broad perspective of proposed QMC scheme for the 
complex solid-state and quantum chemistry problems. 

We consider a fermions system with pair interaction in 
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the most general form and present the partition function 
Z = TrTexp(— S") in the terms of the effective action S: 

S ^ So + W = J Jtlclc'-drdr' + (1) 

Here T is a time-ordering operator, r — {i,s,T} is a 
combination of the discrete index i numbering the single- 
particle states in a lattice, spin index s =t or J, and the 
continuous imaginary-time variable r. Integration over 
dr implies the integral over dr, and the sum over all lat- 
tice states and spin projections: J dr = J^i lo ^.t. 
We borrow the linear-algebra style for sub- and super- 
scripts to make the notation clearer. The creation (cj) 
and annihilation (c'') operators for a fermion in the state 
r are labelled as covariant and contravariant vectors, re- 
spectively. The labelling for coefficients t, w is chosen to 
present all integrands like scalar products of tensors. An 
additional quantity a^, is introduced for the most effec- 
tive splitting of S to the Gaussian part (5*0) and inter- 
action (W). The parameters a^, are to be chosen later 
to to optimize the algorithm and to minimize the sign 
problem. 

We consider £"0 as an unperturbed action and switch to 
the interaction representation. The perturbation-series 
expansion for Z has the following form: 

Z = ET=oIdriJdr[...Jdr',,n,{n,r[,...,r',,) (2) 

O, — 7„tlDlii/^^''^ ■ . ,„''2fc-i''2fc nJ-ir2...r2fc 
ilk — ^0 fei Wrir2 ■■■ Wr2k-ir2k-L^r' r' r' 

' 1' 2---' 2fc 

where Zq is a partition function for the unperturbed sys- 
tem and 

=<T(ct,c'-i-ap).....(c^, c'^--«;-)>. (3) 

2k '1 '1 '2k '2k 

Hereafter the triangle brackets denote the aver- 
age over the unperturbed system, < A >= 
ZQ"^Tr{TAexp(— 5*0)}. Since the action Sq is Gaussian, 
one can apply the Wick theorem and find the expression 
for D in terms of a determinant of 2k x 2k matrix: 

^ir2---r'2k "^r'. V r'.W W 

Here g^, =< Tc^id" > is the the single-particle two-point 
Green function in the QMC notation and Sij is a delta- 
symbol. 

In the following we use an important-sampling Markov 
process in the configuration space, where the points are 
determined by the perturbation order k and the set 
{ri,r[, ...jV^f.}- Suppose for a moment that is al- 
ways positive, and consider a random walk with a prob- 
ability of Z~^Vlk{ri,r'i, ...,r2k,r'2f.) to visit each point. 
Denote the average over this random walk by the over- 
line. Then for example the Green function GJ!, = Z~^ < 
can be expressed as gl,{ri, r[, rl^j^) where 



g^, determines the Green function for a current realiza- 
tion. It is important to note that a Fourier transform 
of g^, with respect to time arguments can be found ana- 
lytically. Therefore the Green function can be calculated 
directly at Matsubara frequencies. Such an approach has 
an advantage over the calculation in r-domain, because 
it automatically takes into account the invariance of the 
initial action in the translations along r-axis. Higher- 
order correlators can be calculated in the same way. More 
detailed description of the algorithm as well as method- 
ological discussion can be found in Ref. [l^ . 

In certain cases proper choice of a can indeed com- 
pletely suppress the sign problem. For example, for 
Hubbard model it is reasonable to choose a-f*,j, — 
asSrr'SijSss' ■ If the Gaussian part of action does not 
rotate spins, than oc Sggi , and the determinant in Q 
is factorized: = 0^0^. For the case of Hubbard model 
with attraction one should choose — ai = a, where 
a is a real. For this choice — g|, and consequently 
Z?l — Di- All terms of are positive in this case, because 
w < 0. 

The choice of a-^ = is useless for a system with re- 
pulsion, because the alternating signs of flk with odd and 
even k appear 0. Similarly to the discrete Hubbard- 
Stratonovich transformations 0], the particle- hole sym- 
metry can be exploited for a half-filled system. One can 
show that a choice aij — 1 — ai = a delivers a condition 
= —Di for this case, thus eliminating the sign prob- 
lem [l^ . Further, for a particular case of an impurity 
problem in the atomic limit = 1 — ai = a with a > 1 
or a < eliminates the sign problem for the repulsive 
interaction at any filling factor p^ . 

Summarizing up these observations, we can write a 
draft recipe of how to choose a. For a physically reason- 
able split of the actionnthe value of a should not be too 
large. Therefore for the diagonal repulsive terms of the 
interaction matrix we propose to use a'\ = 1 — ai = a 
with a slightly above 1. For the attractive interaction, 
and for the off-diagonal matrix elements of w, the choice 
should be q!| = ~ 0.5. Of course, in a general case 

is not positive-defined and one needs to work with its 
absolute value in QMC sampling. In this case an ex- 
ponential fall-off occurs for the large systems or small 
temperature. It is worth while to mention that an above 
choice of parameters a suppress the sign problem for local 
DMFT-like action with diagonal in orbital indices bath 
Green function. 

Now we discuss how to organize a random walk in 
practice. We need to perform a random walk in the 
space of fc; ri, rj^, Tjj,. Two kinds of trial steps are 
necessary: one should try either to increase or to de- 
crease A; by 1, and, respectively, to add or to remove 
the four corresponding operators. A proposition for 
T2k+i,r2k+i , ''2fc+2, ''2fc+2 should bc generated for the "in- 
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FIG. 1: Local Green function of the two-band rotationally in- 
variant model at the Matsubara frequencies. Filled and open 
circles correspond to the static and to the nonlocal in time 
spin-flip, respectively. High-frequency asymptotics is drawn 
with line. Inset shows the distribution function for the per- 
turbation order k. 

cremental" step. The normalized modulus 

\Hr'\w:t:irii:i\ (s) 

\\w\\=JJJJ\w;^'\drdRdr'dR' 

can be used as a probability density for this proposition. 
Then the standard Metropohs acceptance criterion can 
be constructed using the ratio 

j~)ri...r2k+2 

k + l' D'^j-l" ■ ^' 

The " decremental" step can be organized in a same way. 

The most time consuming operation of the algorithm 
is a calculation of the ratio of determinants, defined by 
the Eq. Q). Fast update trick can be used, resulting in 
oc fc^ operations P, 0. Here we estimate k. An aver- 
age value of © determines an acceptance rate for QMC 
sampling. It is reasonable to expect that by the order 
of magnitude this rate is not much less then unity. The 
ratio of determinants times | |w| | can be interpreted as an 
expectation value for \W\. Therefore 

k^W\- (7) 

For the Hubbard lattice of N atoms with an interaction 
constant [/, for instance, \W\ oc P\U\N. In principle, one 
can manipulate with a to minimize \W\. These manipu- 
lations should however preserve the average sign as large 
as possible. 

We apply the proposed continuous time QMC for the 
important problem of super-symmetric two band impu- 
rity model at half-filling |2iJ,|2l|. To our knowledge, this 
is the first successful attempt to take the off-diagonal 
exchange terms of this model into account. These terms 
are important for the realistic study of multi-band Kondo 
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FIG. 2: Imaginary-time Green Function (upper panel) and 
the four-point correlator x (lower panel) for two band model. 
Upper inset shows DOS computed from the Green function. 
Solid and dot lines correspond to the static and to the nonlocal 
in time spin-flip, respectively. Lower inset shows DOS for 5 
band model with the same value of U and J=0.2 

problem, because they are responsible for the local mo- 
ment formation l20|. The interaction in this model has 
the following form 

^(7V(r)-2)(7V(r)-2)-^(S(r).S(r) + L(r).L(r)), (8) 

where N is the operator of total number, S and L are to- 
tal spin and orbital-momentum operators, respectively. 
The interaction is spin- and orbital- rotationally invari- 
ant. The Gaussian part of the action represents the diag- 
onal semicircular density of states with unitary half- 
band width: t{uj) = '^/{oJn + -v/o;^ — 1), where LOn are 
Matsubara frequencies related to imaginary time vari- 
able. We used parameters U = 4, J = 1 at /3 = 4. A mod- 
ification of this model was also studied, where spin-flip 
operators were replaced with the fully non-local in time 
terms. For example, operator cj|.^c°''-'^cj|^c^^'^ was re- 
placed with j3~^ J dr' Cq^^c'^^'^ c[^^,c^'^'^ . Figures present 
the result for the local Green function GH and the four- 
point correlator x(t — t') =< Cq^^c'^^'^c\^^,c^^'^ >. The 
later quantity characterizes the spin-spin correlations and 
would vanish if the exchange is absent. 

Figure 1 shows the Green function at Matsubara fre- 
quencies. The typical number of QMC trials was 2 • 10''. 
Results for the local and non-local in time spin-flip are 
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shown with filled and open circles, respectively. The dis- 
tribution function for the perturbation order k is drawn 
in the inset. For the system studied it appears to be a 
Gaussian-like peak located at A; « 75, with an acordance 
to Eq. ([7|). The estimated error bar in G(w„) is about 
3 • 10~^ for the lowest frequency and becomes smaller as 
the frequency increases. The high-frequency tail obeys 
an asymptotic behavior —lm.(iw + e)~^ with e w 2.9. 

Green function in the time domain was obtained by 
a numerical Fourier-transform from the data for G'(u;„). 
For high harmonics the above-mentioned asymptotic was 
used. Results are presented in the upper panel of Fig. 2. 
The lower panel presents the result for x(t). Thess data 
are obtained similarly, the difference is that x(<^) is de- 
fined at Bose Matsubara frequencies and obeys a l/w^ 
decay. It is interesting to note that Green function is 
almost insensitive to the details of spin-flip retardation. 
Both Green functions are very similar and correspond 
to qualitatively the same density of states (DOS). The 
maximum-entropy guess for DOS is presented in the in- 
set to Fig. 2. On the other hand, switch to the non-local 
in time exchange modifies xi''') dramatically. The local 
in time exchange results in a pronounces peak of x{t) at 
T « 0, whereas the non-local spin-flip results in almost 
time-independent spin-spin correlations. For realistic de- 
scription of Kondo impurities like cobalt atom on metal- 
lic surface it is of crucial importance to use the spin and 
orbital rotationally invariant Coulomb vertex in the non- 
perturbative investigation of electronic structure. The 
proposed continuous time QMG scheme is easily general- 
ized for a general multiband case. As example we shows 
the DOS for five d-orbital model at half-filling for the 
same value of U and J=0.2 in the lower insert of Fig. 2. 

For a final discussion it is suitable to analyze a conver- 
gence of the series (01 . Fermi statistics and a finite size 
of the system insure us that the configurational space 
of the problem is of a finite order. Because the per- 
turbation operator W has a finite norm, its powers 
therefore grow slower than fc!. Consequently, from the 
mathematical point of view there is no doubt that the 
series (|2J always converges. Physically it is important to 
note that this convergence is related both with a choice 
of the type of serial expansion and with the peculiarities 
of the system under study. First of all, series 10) con- 
tains all diagrams, including non-bounded. In the ana- 
lytical diagram-series expansion non-bounded diagrams 
drop out from the calculation and the convergence 
radius for the diagram-series expansion differs from that 
of (O. Further, Fermi statistics is indeed important. An 
analog of Q for Bose field can diverge even for a single- 
atom problem 22|, because in this case one deals with 
an infinite-order Gilbert space. It is important to keep 
this in mind for possible extensions of the algorithm to 
the electron-phonon systems and to the field models, as 
these systems are also characterized by an infinite-order 
phase space. A general time-dependent form of the ac- 



tion (Eq. (0)) allowed us to use renormalization theory 
for the Hubbard-like model: in this case local DMFT 
would be a starting point for lattice calculations in order 
to reduce the effective interaction and minimize the sign 
problem. 

In conclusion, we have developed a fermionic continu- 
ous time quantum Monte Carlo method for general non- 
local in space and time interactions. We demonstrated 
that for a Hubbard-type models the computational time 
for a single trial step scales similarly to that for the 
schemes based on a Stratonovich transformation. An 
important difference occurs however for the non-local 
interactions. Consider, for example, a system with a 
large Hubbard U and much smaller but still important 
Coulomb interatomic interaction. One needs to intro- 
duce N"^ auxiliary fields per time slice instead of N to 
take the long-range forces into account. On the other 
hand, the complexity of the present algorithm should re- 
main almost the same as for the local interactions, be- 
cause \W\ does not change much. This should be useful 
for the realistic cluster DMFT calculations and for the 
applications to quantum chemistry 5] . It is also possible 
to study the interactions retarded in time, particularly 
the super-exchange and the effects related to dissipation. 
This was demonstrated for an important case of the fully 
rotationally invariant two band model and its extension 
with non-local in time spin-flip terms. 
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